Description of Pairing correlation in Many-Body finite systems with density 

functional theory 
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Different steps leading to the new functional for pairing based on natural orbitals and occupancies 
proposed in ref. [D. Lacroix and G. Hupin, arXiv:1003.2860 are carefully analyzed. Properties of 
quasi-particle states projected onto good particle number are first reviewed. These properties are 
used (i) to prove the existence of such a functional (ii) to provide an explicit functional through a 1/N 
expansion starting from the BCS approach (iii) to give a compact form of the functional summing 
up all orders in the expansion. The functional is benchmarked in the case of the picked fence pairing 
Hamiltonian where even and odd systems, using blocking technique are studied, at various particle 
number and coupling strength, with uniform and random single-particle level spacing. In all cases, 
a very good agreement is found with a deviation inferior to 1% compared to the exact energy. 

PACS numbers: 74.78.Na,21.60.Fw,71.15.Mb,74.20.-z 
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I. INTRODUCTION 

Nuclear systems [l], H[ or ultrasmall metallic grains [|[ 
offer the possibility to get insight in finite pairing cor- 
relations of systems with varying particle number. The 
introduction of a simple many-body wave packet ansatz 
more than 50 years ago by Bardeen, Cooper and Schri- 
effer (BCS) 0' was a major breakthrough for the under- 
standing and the description of superconductivity. To 
illustrate the advantages and drawbacks of the BCS the- 
ory, in figure [TJ the condensation energy, i.e. the dif- 
ference between the Hartree-Fock (HF) energy and the 
energy of the system obtained with BCS (dashed line) is 
compared to the exact result (solid line) for the picked 
fence pairing Hamiltonian (for details see section HIT)) 0- 
7]. One of the great advantage of the BCS or Hartree- 
Fock Bogolyubov (HFB) theory is the possibility, under 
the price to conserve particle number only in average, 
to grasp part of the correlation beyond the Hartree-Fock 
level while keeping the theory relatively simple. As can 
be seen from figure [TJ the BCS prediction becomes closer 
to the expected result as the number of particle increases. 
Indeed, the BCS theory is shown to be exact in the ther- 
modynamic limit. Besides these interesting aspects, BCS 
or HFB suffer from a threshold at low coupling. In fact, 
when the coupling strength is much smaller than the av- 
erage level spacing between single-particle states, BCS 
identifies with HF while, in reality, correlations built up 
as soon as the two-body interaction is plugged in. In ad- 
dition, even above the threshold, part of the correlation 
are systematically missed. 

The BCS or HFB theories are nowadays standardly 
used in nuclear physics, for instance, within the Energy 
Density Functional (EDF) approach 8, 9] leading to the 
so-called "Single-Reference" (SR-EDF) or "mean-field" 
level of EDF. These tools already provide a rather good 
reproduction of gross nuclear properties. For instance, 
masses can be estimated with a typical precision of 500- 
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FIG. 1. (Color online) Exact condensation energy (red solid 
line) obtained for the picked fence pairing Hamiltonian as a 
function of the coupling strength for 16 (top) and 8 (bottom) 
particles. In both cases, the BCS (green dash line), the pro- 
jected BCS with a projection made before (open blue circle) 
or after the variation (open violet triangle) are also shown. 
In the right, occupation numbers of the different theories are 
plotted for g/Ae = 0.82. 



600 keV. Figure [TJ however clearly points out that there is 
room for improving the BCS approach in finite size sys- 
tems. In particular, part of the discrepancy stems from 
the use of a trial wave-function that is not an eigenstate 
of the particle number operator N. Starting from the 
BCS wave-packet, a new state with good particle number 
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can be obtained using projection operator technique [lOj. 
Within EDF, similarly to the restoration of angular mo- 
mentum or calculation including dynamical fluctuations 
associated to configuration mixing, projection onto good 
particle number enters into the class of Multi-Reference 
EDF (MR- EDF). If the projection is made prior to the 
variation (Variation After Projection [VAP]), the varia- 
tional state directly becomes an eigenstate of N. Illustra- 
tion of VAP condensation energy (open circles) is given 
in figure [T] (see for instance flljl). Such an approach pro- 
vides a very accurate description of pairing correlation at 
all coupling strengths and completely removes the BCS 
threshold problem. VAP still remains rather involved nu- 
merically and a less efficient but simpler approach con- 
sists in projecting the state after the variation, the so- 
called Projection After Variation [PAV] (open triangles 
in figure [I}. Projection technique is becoming a popular 
tool in nuclear structure. However, recent studies have 
shown that projection aiming at restoring broken symme- 
tries and/or more generally configuration mixing should 
be handled with care when combined with density func- 
tional theory [12|, |l3(, due to the possible appearance of 
jumps and/or divergences in the energy surface. These 
difficulties have been carefully analyzed in refs. (T^4l6j 
and have been related to the self-interaction and self- 
pairing problem. By comparing theories starting from 
an Hamiltonian and an energy functional, a correction 
to the pathologies was proposed such that systematic 
calculation along the nuclear chart is now within reach. 
These studies have clearly pointed out that specific as- 
pects might appear due to the use of functional theories 
(see also [HI, El]) wnen MR-EDF is used. 



The EDF framework provides a unified framework not 
only for nuclear structure but also for nuclear dynamics 
and thermodynamics. While MR-EDF is a suitable tool 
for the former, due to its complexity, it can hardly be 
used in the latter cases. The goal of the present work 
is to discuss a new approach to treat pairing where the 
projection effect is directly incorporated into the func- 
tional through specific dependencies on natural orbital 
occupancies. Such an approach, directly written in the 
functional framework, avoids some ambiguities encoun- 
tered in current EDF and is expected to greatly simplify 
both PAV and to be easily adapted to non-equilibrium 
evolution of finite temperature studies. Main aspects of 
the new functional theory have already been summarized 
in ref. [l9| . Here, we present a complete discussion of the 
different steps leading to the functional. Below, we first 
discuss the interest of using natural orbital based func- 
tionals. Then, mathematical properties of Projected BCS 
states that are used to propose the functional, are given. 
Finally, the new functional is applied to a specific pairing 
Hamiltonian either with equidistant or non-equidistant 
level spacing and benchmarked for any coupling strength 
and particle number. 



A. Functionals based on natural orbitals and 
occupancies 

The possibility to replace a many-body problem by 
a functional of the density matrix has been first pro- 
posed by Gilbert in ref. [2(| and is named Density 
Matrix Functional Theory (DMFT) or Reduced DMFT 
(RDMFT). The Gilbert theorem is a generalization of 
the Hohenberg-Kohn theorem [2l[ where the variational 
quantity, i.e. the local density p(r, r) is replaced by the 
full one-body density matrix (OBDM) p(r, r'). Most of- 
ten, the OBDM is first written in the natural or canon- 
ical basis as 7 = JV |y>j)nj(y>i|. Here n, and 
denote occupation numbers and natural orbitals respec- 
tively. Then, the initial many-body problem is replaced 
by the minimization of an energy functional 

?[{<Pi}, K>] = £[{<Pih {»*}] - v{Tr(Np) - N} 

-^2 x ij(i i Pi\ l Pj) (1) 

where the variation is made with respect to both single 
particle states <p*(v) and occupation numbers. The set 
of Lagrange multipliers fi and {Ay} are introduced to in- 
sure particle number conservation and orthogonality of 
the single-particle states. RDMFT has several advan- 
tages compared to standard Density Functional Theory 
(DFT). For instance, while Kohn-Sham single-particle 
states used to construct the local density are not ex- 
pected to have physical meaning, the non-local density 
7 should match the exact one at the minimum. Accord- 
ingly, associated single-particle states and occupations 
identify with the one of the exact many-body state. This 
is an important aspect of this theory. Indeed, DFT can 
only provide information on the energy. In RDMFT, not 
only the energy can be estimated but also any one-body 
operators. Similarly to Density Functional Theory, the 
main challenge is to find accurate functionals. 

Another interesting feature of this theory is its abil- 
ity to describe aspects that arc not adequately obtained 
at the DFT level, like reactions, atomization energy or 
the dissociation of small molecules. All these phenomena 
have their counterpart in nuclear physics. Nowadays, a 
sizeable effort is made to provide new accurate RDMFT 
functionals and benchmark them on finite and infinite 
systems (see for instance [22j and refs. therein). 

In this article, we focus on pairing. Let us first re- 
mark that current SR-EDF that account for pairing al- 
ready share many aspects with RDMFT. Most nuclear 
SR-EDF used nowadays start from a functional that can 
be written as 

£ SR [ P ,C] = £ P + £ PP + £ C 

ij ijkl 

+ \Y J vg ij C ijM . (2) 

ijkl 
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where v pp and v c denote effective two-body kernels re- 
spectively in the particle-hole and correlation channels. 
C± t 2 denotes the irreducible two-body correlation matrix 
defined as the difference between the two-body density 
and the antisymmetric product of one -body density ma- 
trix (see for instance |23j). To treat pairing correlations, 
a quasi-particle trial state, \4>qp), is considered, then the 
correlation matrix elements can be written in terms of 
the anomalous density k as Cij.ki = KijKki [H, (3- I n 
the natural orbital basis, the quasi-particle state can be 
expressed in a BCS form 



\M=Y[(l + Xl alal) |0>, 



(3) 



where |0) corresponds to the particle vacuum while 
{aJ,aJ} correspond to doubly degenerated canonical 
states {^Pii^pi} with occupation probability 2rii. The 
Xi coefficients are connected to the occupation numbers 
through 



UA 2 = 



(1-rnY 

Accordingly, pairing energy reduces to: 



(4) 



probabilities of the new trial state \N) do not appear di- 
rectly [24]. Nevertheless, the occupation numbers of \N) 
can be estimated numerically. Illustration of occupation 
probabilities of the projected states are compared to the 
exact ones for the picket fence pairing Hamiltonian in 
figure [TJ Both VAP and PAV results as well as BCS case 
are displayed. It is first interesting to mention that, while 
the energy is improved in the PAV case, single-particle 
occupation numbers deviate more from the exact solu- 
tion than the original BCS case. This is something to 
worry about since, when PAV is performed in EDF, ex- 
pectation values of one-body operators are estimated. In 
opposite, in the VAP method, occupation probabilities 
perfectly match the exact case for all particle number 
and pairing coupling strength. Therefore, we see that 
the use of projected state before the variation leads to a 
very good reproduction of both the ground state energy 
and the single-particle occupation numbers. 

Having this in mind, in the following, we use the prop- 
erties of PBCS state to provide a new functional for pair- 
ing directly based on the occupation numbers of the pro- 
jected state. The state (JSJ) is used as a starting point 
where it is implicitly assumed that the orbitals are writ- 
ten in their canonical basis, namely the one which ex- 
hibits an explicit time reversal symmetry. In that case, 
the energy © reduces tcQ 



Noting in addition that both £ p and £ pp can directly be 
written as a functional of and (ft through their depen- 
dence on the one-body density, we see that current SR- 
EDF can indeed be interpreted as a mapping between the 
initial problem into a functional theory of ({v?i}> {fit}), 



ij 



v pp o N o N 

ij 



£[ P ,C] 



£[{<Pi},{ n i}] 



(5) 



provided that the functional is written in the canoni- 
cal basis. EDF based on quasi-particle states have the 
shortcomings discussed in the introduction when com- 
bined with projection onto good particle number within 
MR-EDF. This is nowadays used in nuclear structure 
study. If the projection is made prior to the varia- 
tion, such a projection is equivalent to consider a new 
trial wave-function, called hereafter generically Projected 
BCS (PBCS) state, of the form: 



TV 



\N) = P n \4> Q p) cx 



xia H 



|0>, (6) 



where P N is the projector on particle number N (see for 
instance [ltj [HI, EB])- In the following, we will use the 
short notation = J2i x i"l with b\ = a\at. 

When projection is made in EDF, the associated func- 
tional becomes much more complex to minimize. It 
should however be noted that the functional is generally 
written in terms of the normal and anomalous density 
of the original quasi-particle state from which the pro- 
jected state is constructed. Therefore, the occupation 



were p 



N 



and C>! .-. now stand for the occupation 
and correlations associated with the projected state, i.e. 



(N\N) 



C 



JV 



(N\b\b 3 \N) 
(N\N) 



6 l] n i rij . 



(8) 



Here, we have used the compact notation Cjj 



",33 



In the following, we will omit the N label to shorten no- 
tations keeping in mind that these quantities refer to the 
projected state. In order to do the mapping ([5]), we are 
left with the challenge consisting in expressing the corre- 
lation Cij as a functional of n, as it can be easily done 
in the BCS or HFB case. But in the present work, we 
aim at accounting for the particle number conservation 
directly in the functional. 



1 Note that, correlation matrix elements should also appear in the 
particle-hole channel. Since the aim of the present article is to 
focus on pairing channel and since these components cancel out 
exactly in the example presented below, they are omitted here. 
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II. CONSTRUCTION OF FUNCTIONALS FOR 
PAIRING FROM A PBCS STATE 

Here, some properties of projected states are first high- 
lighted. These properties are then used as a guidance to 
construct the functional. Over years, interesting features 
of matrix elements entering in Eq. (jHJ have been derived. 
Some can eventually be deduced using the fact that the 
BCS state plays the role of PBCS state generating func- 
tion [25l l2q | and can be used, for instance, to minimize 
the energy directly written as a functional of the {x{\ pa- 
rameters [13, HI] . Proofs of some of the properties that 
are used below are first given. 



A. Definition of a class of operators, states and 
overlaps 

First, we start with a strategy similar to ref. [2j|. A 
set of pair creation operators that omit one, two,... pairs 
of single-particle states is first introduced: 

= rt - Xl b\ 



r+( 



r\ij) = - xib\ - Xj b], 



(9) 



where indices i, j refer to the removed pairs. In the fol- 
lowing, will denote the size of the single-particle Hilbert 
space From these operators, a corresponding set of states 
with a given particle number is defined: 



\K) 



c K (rt) 



K 



\K:i) = c K {tf(i)) K \-) 



\K:i,j) = c K (rt(z,j)) K |-) 



(10) 



with K < N while ck is taken by convention equal to 
(Kiy 1 / 2 . Note that, the state introduced in Eq. © 
corresponds to the special situation where K = N and 
no pair has been been removed. From these states, we 
define a set of coefficients from the overlaps: 



I K K\ (K\K) 

I K (i) = Kl (K : i\K : i) 
I K {i,j) = K\ (K:i,j\K:i,j) 



(11) 



Using the fact that (frj) 2 = 0, due to the fermionic nature 
of the particles, the different operators verify: 

(ft (h, ■ ■ ■ ,i K )) K = (r f («!, • • • ,i K ,j)) K 



+ K Xj (^(h,--- ,i K ,j)f \(12) 

This property leads to specific relationships between the 
states defined above and their overlaps. For instance: 



l K 



= I K (i)+K\xi\ 2 I K -i(i), 



I K (i) = lK{i,j) + K\x \ 2 I K -i{i,j)- 



(13) 



These recurrence relations have been recently used to 
solve numerically VAP (2?| and will be at the heart of 
the present work to design a new functional for pairing. 



B. Energy as an explicit functional of {xt} 

Since the PBCS state is written as a functional of the 
parameter set {xi}, expectation values of any operators 
can a priori be expressed as a functional of this set. Here, 
an illustration is given for the occupation probabilities 
and correlation matrix elements. 

Using the states defined in Eq. (fT0|) . expectation values 
of operators entering into Eq. (jSJ can be expressed as 

(N\a\a l \N) = \x l \ 2 (N-l:i\N-l:i) 

(N\blb 3 \N) = x* Xj (N -l:i\N-l:j). (14) 

We then deduce that both occupation numbers and cor- 
relation components can be expressed in terms of ratios 
between the different coefficients introduced in Eqs. (fTTT) : 



N\x, 



2 Jjy_i(i) 
In 



C« = Nx;x s Iff -^ j) for (i^j) 



(15) 



while for i = j, Cu = n;(l — m). Overlaps entering in rii 
and Cij can be directly expressed as a functional of {xi}. 
Indeed, a direct development of (T*) K in © gives: 



\N) = c K 



Xj-, ' ' ' Xj^j b-_ b: 



= K\c K 



E 

i<— <ijr<n 



x%x • • • Xi N b it ■ ■ ■ b iN j 



where . >, is used to insist on the fact that the 

summation is made only for indices different from each 
others. From this expression, it is straightforward to see 
that 



E 

(ii,— ,ifc) 



(16) 



In a similar way, the following expressions can be de- 
duced: 



f Ik® 



/-((ii,— ,i K 



)jti \ x ii I 



lK(i,j) - E(t 



>*K)#(i,j) 



Note that, these expressions also suggest additional re- 
currence relation between the overlap: 



Ik = E^\ 2 lK(i) 
Ik{i) = T,&i\xj\ 2 lK(iJ) 



(17) 
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For completeness, additional properties are given in ap- 
pendix [A3 Reporting above expressions into ([T5"j) . both 
Ui and Cij , and consequently the energy, take the form of 
an explicit functional of {xi}. This functional turns out 
to be too complex for a direct practical use unless one 
can take advantage of the different recurrence relation to 
estimate the desired quantities (27j . 



and {x^ is of particular interest for the present discus- 
sion. Indeed, given a set of occupation numbers n i; one 
could a priori deduce the values of the Xi through these 
secular equations. This shows that these parameters are 
implicit functional of the occupation probabilities (see 
also discussion in section HTT|) . 



C. Energy as an implicit functional of {n^} 

The possibility to write the energy as a functional of 
natural orbitals and occupation probabilities is far from 
being trivial. Strictly speaking, Gilbert theorem [2(| 
holds for systems bound by an external potential. It 
could however be extended to self-bound systems with 
the introduction of Legendre multiplier technique [3(| ■ In 
practice, such a technique is useful when the energy can 
first be written as a functional of the single-particle en- 
ergies throu gh s ome preliminary approximations (see for 
instance (3ll |32|). In general, the existence of occupation 
number functional as well as its form is not straightfor- 
ward. Here, we give a proof of principle that the energy 
estimated with a PBCS trial wave can indeed be written 
as such a functional. Since all quantities can be written 
as a functional of the {xi}, it is sufficient to prove that 
these parameters can in turn be put as a function of the 
{rii} set. 

Starting from the expression of rij and taking advan- 
tage of (p"7j) . we first obtain: 



n i =N'£\x i \*\x j \* Is -- 1 



In 



(18) 



Then, using the following recurrence relations 

I N -i(i) = In-i(i,]) + (N -l)\x 3 \ 2 I N _ 2 (i,]) 

I N -i(j) = I N -i(i,]) + (N -l)\ Xl \ 2 I N _ 2 (iJ), 
which are valid for any i ^ j, we see that: 



lN-l(i,j) 

lN-2(i,j) = 
from which we deduce 

n t (N 



\xj\ 2 I N -i(j) - \x l \ 2 I N - 1 {i) 



LN-l 



(0 



IN -l 



(j) 



(19) 



N - 1 



\xA 2 rii 



Eventually, it can be transformed as: 



\xi\ 2 nj 



W(l -«<)=X>i-ni) r^ji- 



(20) 



(21) 



This expression holds for any single- particle state i. This 
set of coupled equations between occupation numbers 



D. Energy as an explicit functional of {n, } 

In this section, we discuss the main objective of the 
present work, i.e. to provide an explicit functional of the 
occupation probabilities. The strategy that is followed 
here is to use the BCS case as a guidance (see Appendix 
iBj). In that case, there is a direct and simple relation 
between \xi\ 2 and rn already given in Eq. (HJ). Let us 
first see how this relation can be generalized in the PBCS 
case. 

Using the first equation of (fTB")) for K — N and report- 
ing in the denominator appearing in ni, leads to 



(22) 



where we have introduced the notation ajv(i) = 
^iv(«) / {NIff-i{i))- This expression can easily be inverted 
and compared to (gj. In the PBCS case, we have: 



1 



a N (i). 



(23) 



Therefore, we see that the BCS limit is recovered if 
ajv(i) = 1 and that all the physics beyond the ordinary 
BCS or HFB theories is contained in its deviation from 
one. This could also be seen by expressing the correla- 
tion in terms of and ccjv(z). Reporting Eq. (|19[) into 
([15]), leads to 



Cij — 



ni{l-m) for (i=j), 



Fil 2 - Fil 



(24) 



Taking advantage of (|2"3"|) and using the short-hand nota- 
tion a, = ajq{i) 1 finally gives (for i ^ j) 



= \ rii(l - 7ij)rij(l - rij)aia 



J V 3 J £ J 

■rii — n-, 



In the limit Ui 



rii(l - nj)^ - rij{l - ni)otj 

1. the BCS functional Cj. 



(25) 



y/rii(l — rii)nj(l — rij) is recovered. More generally, it 
is shown that any of the following quantities, defined 
through: 



a K (ii 



1 Ik{h, 



K I K -i(ii,--- ,i K ) 



(26) 



identify with 1 in the BCS limit (see appendix |B|) . 
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1. 1/N expansion beyond the BCS theory 

Since the BCS theory identifies to PBCS in the large 
N limit, it is reasonable to seek for a correction to 
oti((ii,'" = 1 written as a 1/N expansion. Such 
an expansion can be obtained thanks to the relation: 

otK(h, - ■ ■ , ik) = 

— V l a I 2 ajt - 1 (' 1 '-'' y -J') (27) 
K ■ , %\ 2 + <*K-i(ii,--- ,i K ,j)' 1 j 

connecting % and (XK-i terms. This expression can 
be derived using (fT3| and (|T7|) . Due to the presence of 
a 1 / K prefactor in this relation, any correction of order 
1/(K—1) in a.K-1 will appear as as an order 1/K(K — 1) 



in ax- As an illustration, assuming that otN—i(i,j) ~ 1 
as in BCS, leads to: 

= I(iV-n i ) = l-ir ii! (28) 

that appears as the first order correction in (1/N) to the 
BCS case. Similarly, we can obtain: 

ajv-i(i,j) ^ jV _ (iV-rtj -nj) 

aN-2(i,j) - (iV - Tti - nj -rifc) 



Higher order corrections in ajv(«) can be obtained by including more and more terms in the expansion of all aj<- 
(with K < N). This technique has been used in [l!| to get the expansion: 



11 1 4. 

a N (i) =l--m+ N{N _ l) E n ) t 1 - ( n > + n j)l + N (N-l)(N-2) fc S ^ [2 ~ (nj + ^ + nfc)] 



(29) 



which corresponds to ccjv(i) written as an explicit func- 
tional of the occupation numbers. Note that, additional 
terms tested numerically as negligible and appearing at 
the second (or higher) order approximation are omitted 
here. This functional can then be injected into (l23t lead- 
ing to an explicit functional of Xi in terms of the rii . Ac- 
cordingly, expectation values of any operators becomes 
also a functional of the projected state occupation num- 
bers. 

This approximation has been tested numerically in ref. 
[l9| and has shown a rapid convergence in the strong 
coupling limit. However, for small coupling (HF limit), 
a slow convergence was found. Indeed, assuming that 
rii — > 1 for the N pairs, we deduce, for one of the occupied 
state: 

£„2[l_(„,.+ nj .)]->-(JV-l) 

J2 n K t 2 - ( n < + u j + n k)] -> -( N - !)( N - 2 ) 



Therefore, in this limit, all contributions to any order 
will participate to the same extend and sum-up to give 
&N(i) = (1 — n i) leading finally to a correlation given by: 



y/rnrij 



(30) 



This form, which has been proposed using a completely 
different strategy in electronic system (33|, will never be 



properly described by the BCS functional. From the dis- 
cussion above, difficulties in the application of the present 
functional might also be anticipated. Indeed, since all 
terms in the expansion should be kept, the functional 
becomes rather complicated and its application might be- 
come rapidly intractable. 



2. Re-summation of the 1/N expansion and simplified 
functional 

The price to pay to correctly describes the weak cou- 
pling limit is to keep all orders in the expansion pre- 
sented above. This basically shows that the 1/N expan- 
sion approach starting from the BCS approximation is 
not appropriate in that case. To overcome this difficulty 
a simplified functional can be found using the following 
approximation in Eq. (|29[) . 



X(N-l) £ ^ AT2 £' 



t 1 



N(N -1)(N -2) ^ N 3 

while keeping all terms in this expansion. This approx- 
imation leads to a simple linear dependence of the oti 
coefficient with respect to the occupation numbers rif. 



a.% = a - airii, 



(31) 



7 



where ao and a\ are given by the expressions: 
1 



III. APPLICATION 



01 = — (1 + s 2 + si + 

1 1-4 

N 1 - s 2 



(32) 



and 



«- = :i + (l + 2s 2 + ■■■ + (N — l)s»- 2 ) 



N 

1 + (S 2 - S 3 ) 



9a i 
ds 2 ' 



(33) 



and where the moments s p = — ^](^i) P have been used. 

Reporting expression pip in correlation matrix elements 
Eq. (|25[) gives the simple form (for i =/= j) 



\Jrii(l- n.j)n 3 (l - n 3 ) 

;< \/ ( Q o ~ Q i n ») ( Q o " airij) 
ao — ai (n, + n,- - mrij) 
C(m, rij) 



(34) 



The functional (|3"4"|) together with (|32j|33p represent the 
main result of this article. We can already anticipate 
some advantages of this functional (i) In the Hartree-Fock 
limit s p = 1 for all p > 1. Accordingly, a = ai = 1 and 
we recover the HF functional quoted above, i.e. Cy = 
y/riiTij. (ii) The BCS limit is also easily identified in ([34]) 
by taking the limit ao = 1 and a\ = 0. The net result 
of our approach is that the energy introduced in Eq. ([7]) 
that was originally written as a functional of the density 
and correlations in the projected state becomes now a 
functional of the one-body density matrix components 
only. In practice, such a functional approach should be 
solved by minimizing (|T|) where the energy now reads: 



pp 

iijj "'" J 



i ij 



First applications of this functional can be found in ref. 
[l9j illustrating the predicting power of the functional 
for energies and occupations probabilities. In numerical 
implementation, sequential quadratic programming leads 
to very good convergence at any coupling and/or large 
particle number. Below, the new functional is further 
illustrated and benchmarked. 



We consider here a system of A particles interacting 
through the pairing Hamiltonian of the form [H-Q 



H 



i{a\ai + a]a-A 



t t 

i,3 



(35) 



where i denotes the time- reversed state of i, both asso- 
ciated with single-particle energy £j. The total single- 
particle Hilbert space size is assumed to be 17 = 2A. 
This Hamiltonian can be solved exactly numerically by 
making use of the so-called Richardson equations. First 
test of the functional have been made with this model 
Hamiltonian for an even particle number A = 2N where 
N denotes the number of pairs and for equidistant single- 
particle levels, the so-called "picked fence" Hamiltonian. 
Here, we will further illustrate some of the aspects of the 
new functional in that case, and extend the application 
to even systems (^4 = 2N + 1) and/or non-cquidistant 
levels. 



A. Illustration in the picked fence Hamiltonian 

Here, we first consider the special case of equidistant 
single-particle levels with a level spacing denoted by Ae. 
First, we remind that the strategy to design a functional 
going beyond the BCS one has been made in three steps: 
(i) The parameters {x{\ have been first shown to be im- 
plicit functional of the {rii} through the existence of a 
set of secular equation (Eq. (f2"Tj) ). (ii) Starting from the 
BCS prescription, systematic l/N corrections have been 
proposed to get a new functional (Eq. ([29]) ) (iii) Sum- 
ming all orders, a simplified functional is then introduced 
(Eq. pTl) ). The step (ii) has been shown to be inade- 
quate [19(, especially in the weak coupling limit. Before 
discussing (iii), the existence of secular equations as well 
as the uniqueness of the relation between {x{\ and {rii} 
sets of variational parameters is analyzed. 



1. Existence and uniqueness of a functional of rii 

Eq. (f2"Tj) is proving the existence of a functional of the 
occupation number, at least an implicit one. A graphical 
illustration of Eq. (f2Tj) at the PBCS energy minimum is 
given in figure O The recurrence method of ref. [27| has 
been used to obtained the PBCS solution in this figure. 
The solid curve corresponds to the right hand side of Eq. 
(|2Tj) divided by (1 — rii) as a function of \x5\ 2 keeping 
other Xi fixed as well as the n^. Horizontal lines corre- 
spond to the \xi\ 2 value at the minimum. The horizontal 



dotted line corresponds to N — 8. Equation (|2ip is ful- 
filled when the dotted line crosses the solid line, which is 
indeed the case for the value of x§ minimizing the func- 
tional. Calculations are done for 16 particles and a pair- 
ing constant g/Ae = 0.22. An open square has been 
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added to underline the physical solution. The main in- 
terest of Eq. (f2Tl) is to prove that the PBCS energy can 
be indeed put, at least implicitly, as a functional of the 
occupation numbers. 



12 




0.5 



1.5 



F5 



FIG. 2. (Color online) Graphical illustration of Eq. (l2Tj) 
for the pairing Hamiltonian with 16 particles (N — 8) and 
g/Ae = 0.22. The solid curve corresponds to the right hand 
side of Eq. (|21|l divided by (1 — ni) as a function of | | 2 
keeping other Xi fixed as well as the n». Vertical dashed lines 
correspond to the \xt\ 2 value at the minimum. The horizontal 
dotted line corresponds to N = 8. Equation (|21[) is fulfilled 
when the dotted line crosses the solid line, which is indeed the 
case for the value of x$ minimizing the functional. An open 
square has been added to underline the physical solution. 

This figure also illustrates that the solid line and the 
horizontal dotted lines cross each other several times and 
one may worry about the uniqueness relationship be- 
tween the {xi} and the {n^}. It should however be kept 
in mind that the crossing highlighted by the open square 
is the only point where Eq. (|21[) is fulfilled together with 
the secular equations for other \xi\ 2 . Indeed, starting 
from the recurrence relation and the expression of the 
occupation numbers (|15|) gives 

n 3 = T~\ x j\ 2 (lN-l(hj) + \Xi?lN-l{i,j)) ■ 
J-N 

Considering two states i and j, we then have 

l^l 2 ) I N -i(i,j). 



In 



(36) 



Since lN-i{i,j)/lN > 0, if rij > rii then \x^ > \xi\ 2 . 
This proves for instance that only crossing points in be- 
tween 1 24 1 2 and \xq j 2 might fulfill the secular equation. 
A careful look at figure [2]) shows however two crossing 
points in this region. Let us assume that two solutions 
| £5 1 2 = a and | a^s | 2 = b might exist and fulfill the secular 
equation while keeping all other Xi fixed. Using above 
equation, it is possible to prove that necessarily a = b 
which finally proves that only one of the crossing is phys- 
ical. 



2. Application of the new functional for equidistant level 
spacing 

We first consider the case of even systems with dou- 
bly degenerated equidistant levels. In the following, the 
condensation energy, denoted by fcond, defined as 



£cjond — £hF — £, 



(37) 



where £hf = %J2i>o £i ~ 9^ ls ^ ne Hartree-Fock (HF) 
energy while £ denotes the energy of the considered the- 
ory, fcjond quantifies the predicting power of different 
approximations. An illustration of the evolution of this 
quantity as a function of the coupling strength has al- 
ready been given in the introduction, Fig. (TTJ). In fig- 
ure [HI we see that the proposed functional is almost on 
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FIG. 3. (Color online) Evolution of the condensation energy 
for the exact (red solid line), BCS (green dash line) and new 
functional (blue filled circles) obtained for the picked fence 
pairing Hamiltonian as a function of the coupling strength for 
16 (top) and 8 (bottom) particles. In the right, occupation 
numbers of the different theories are plotted for g/Ae — 0.82 
(top) and g/Ae = 0.22 (bottom). 



top of the exact result (and the exact VAP calculation). 
A slight difference is observed in the intermediate cou- 
pling regime. Similarly, occupation numbers perfectly 
match the exact ones in the strong coupling regime and 
slightly differ from them below the BCS threshold. In 
this regime, while BCS identifies with HF, here, occupa- 
tion probabilities different from 1 and are obtained as 
soon as the interaction is switched on. 



9 



3. Critical discussion of the linear approximation (Eq. 

Wt) 

Figure [4] shows the accuracy of the present approxima- 
tion in the model case of a constant two-body interaction 
g. In this figure, the approximate a.; L for different coupling 
strength g/Ae = 0.32 (filled circles), 0.64 (crosses) and 
0.96 (open circles) are compared to the exact ones, (re- 
spectively dashed, dotted and solid lines) as a function of 
either the orbital probabilities (left) or single-particle en- 
ergies (right) at the minimum of energy. The dependency 
of a,i obtained in the PBCS case for small coupling also 
shows that a simple linear approximation cannot fully 
grasp the physics of weak coupling. Following the same 
strategy as above, quadratic or cubic corrections might 
eventually be obtained. However, this will add complex- 
ity to the functional while the energy is already rather 
well reproduced. 



1.5 



0.5 





1 1 1 _ L. —1 . 








1 


\ 




****->< ^ 












\ 




\ • 


A 






i.i.i. 


i.i.i 



0.4 0.8 5 10 15 

Ui e i/A£ 



4£S 



9 








6 








3 


7/ 

/' 

/' A 
PI 
/ ' 

y i 


3 



1 ■ 

/ 

ni lT_ / A 

■ . i 







_ CC / 




0.2 0.4 " 




0.3 


0.6 


0.9 1.2 



7a. 



FIG. 5. (Color online) Evolution of the one-body entropy for 
different theories as a function of the coupling strength for 
16 particles. The BCS (green dashed line), PAV (open violet 
triangles) and PBCS-functional (filled blue circles) ansatz are 
compared to the exact result. The inset magnifies the low 
g/Ae vicinity. 



threshold, the new functional is in close agreement with 
it. 



5. Simplified, functional for the strong coupling regime 



FIG. 4. (Color online) Evolution of the coefficients cm as 
a function of m (left) or Sj (right) at the minimum of en- 
ergy. The different curves correspond to the PBCS result for 
g/Ae = 0.32 (dashed line), 0.64 (dotted line) and 0.96 (solid 
line). The corresponding results obtained with the linear ap- 
proximation (Eq. (J3TJ) are displayed by filled circles, crosses 
and open circles respectively. 



4- Systematic analysis of occupation numbers 

In figure [31 illustrations of occupation numbers obtain 
in different theories are shown for specific couplings. In a 
DMFT framework, not only the energy should match the 
exact energy at the minimum but also the deduced one- 
body density matrix and a fortiori occupation numbers 
should also be identical to the exact one. To systemati- 
cally compare the gain in predicting single-particle occu- 
pation numbers in the new functional, we have plotted 
the one-body entropy 

S[ n i] = - ^2i n i lo g(™0 + (! - Th) iogl 1 - n i)} ( 38 ) 

i 

in figure [SJ While PAV and BCS are unable to reproduce 
the exact result, especially below or in the vicinity of the 



One important issue from the practical point of view 
is the possibility to further simplify in some regime. In 
particular, in the strong coupling regime, we have seen 
that 1/N perturbation starting from BCS rapidly con- 
verge to the exact solution. Truncation at second order 
of eq. dHJ) already gives a very good result [lj|. It is 
therefore legitimate to question whether a simpler form 
for can be found in this regime. Close to BCS, we 
expect ao — > 1 and a\ — > which plaid in favor of an 
expansion in orders of (ai/ag). For instance, 



Cij = yjmnj(l - Hi)(l - n 3 ) 

In Figure [6[ leading order (LO) [top] and next to next 
to leading order (N 2 LO) [bottom] are compared as a func- 
tion of the pairing strength for a typical number of par- 
ticle A = 16. from and 1 to the exact distribution. It 
can be inferred that the full functional solution can only 
be recovered at low coupling strength when all terms of 
the expansion are taken into account which agrees with 
the previous discussion leading to resummation. 
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FIG. 6. (Color online) Comparison between leading order 
(top) and next to next to leading order (bottom) of the new 
functional (blue filled circles) for 16 particles. Evolution of 
the condensation energy for the exact (red solid line) and 
BCS (green dash line) are shown as references. 



6. Application to odd systems 

Similarly to the BCS framework, the energy of systems 
with an odd number of particles can be obtained by using 
blocking techniques. In the PBCS case, this is equivalent 
to consider a modified trial wave function given by 

|2tf+l)oc(at +4) (r+( a ))V> (39) 

which do preserve the time-reversal symmetry of the so- 
lution. Here, {a, a} correspond to the blocked pair, and 
identify with the last occupied levels in the Hartree-Fock 
limit. The particle number conservation implies that oc- 
cupation of the blocked states are kept fixed and equal 
to rib = ni = 0.5, which is nothing but the filling ap- 
proximation for doubly degenerated states. As an illus- 
tration of the odd-even effect, we define the average gap 
A through the relation: 

A = £c (40) 

T,ijib V^(l-ni) 

This quantity identifies up to a factor l/g with the stan- 
dard gap in the BCS limit. In figure [7J the evolution 
of A/A as a function of particle number A is presented 



for different values of the coupling strength in the ex- 
act (solid line), BCS (dashed line) and new functional 
(filled circles) cases. On the left, odd particle numbers 
are shown as compared to even ones (right), so to dis- 
tinguish odd-even effect. This figure shows that the new 
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FIG. 7. (Color online) Evolution of A/A as a function of 
particle number A for even (left) and odd (right) systems. 
From top to bottom, the three different coupling constants 
g/Ae = 0.66, 0.44 and 0.224 are shown. In each case, the BCS 
(green dashed line) and PBCS (blue filled circles) functional 
theories are compared to the exact calculation (red solid line). 
Note that g/Ae = 0.224 is below the BCS threshold for some 
values A which leads to an equivalent threshold in the quan- 
tity A/A. In the insets, the standard deviations to the exact 
calculation renormalized to 1 are compared for the different 
functionals. 

functional predicts well A/A for both even and odd num- 
ber of particles. Deviations at low coupling strength of 
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the PBCS from the exact case stem from the small dis- 
crepancies in the occupation numbers between those ob- 
tained in the functional formulation and the exact ones. 
The insets of Fig. [7] show the standard deviations from 
the exact calculation normalized to unity for the different 
functionals. It is worth mentioning that the same accu- 
racy is observed for both even and odd systems in the 
case of the PBCS functional, this is in contrast with the 
BCS calculations. In the following discussion, the effect 
of particle number is further investigated. 

7. Accuracy of the functional with respect of particle 
number 

It is known from [34| that the PBCS state exhibits 
slight deviations from the exact solution for medium 
number of particles. Since our approach is based on a 
PBCS trial state, we do expect a similar behavior. To 
systematically address the quality of the PBCS func- 
tional with respect to both the number of nucleons and 
the coupling strength, the condensation energy for odd 
and even systems are displayed in figure [8] as a function 
of [H d/A where 

d/A=j^xh(l/g) 

for g/Ae = 0.224 (top) and g/Ae = 0.44 (bottom). In 
this figure, particle number ranging from A = 8 (large 
d) to A = 360 (small d) have been used. This figure 
illustrates the improvement of the new functional com- 
pared to BCS. It also clearly shows, that some deviations 
from the exact results persist in the new functional. It 
should however be kept in mind that the observed devi- 
ations correspond to less than 1% of errors in the total 
energy. This is illustrated in the insets of figure [8j where 
the relative error defined through 



'-■exact 

where fexact is the exact energy, is displayed as a function 
of d/A. As expected, error tends to zero in all cases as 
A increases (tf — » 0). For intermediate to high coupling 
(Fig. bottom), a good agreement between the PBCS 
based functional and the exact solution is obtained, while 
at lower coupling strength some deviations appear. This 
results both from the approximation scheme used to de- 
sign the functional (linear approximation for the aj, see 
[19| ) and from the accuracy of PBCS theory itself as an 
approximation of the exact trial wave-function. It should 
indeed be kept in mind that the present functional is 
entirely based on the PBCS theory which already devi- 
ates from the exact solution (see for instance |35|). As 
a consequence, it could only lead to results which are 
at most equivalent to the PBCS approximation. From 
the comparison between Fig. [5] (top) and ref. [36j], it 
can be inferred that deviations at low g stem from (i) 
The deviation of PBCS result from the exact solution as 
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FIG. 8. (Color online) Condensation energy predicted by the 
BCS functional (green dash line) and the PBCS functional 
(blue filled circle) compared to the exact solution (red solid 
line) for varying d/A and for g/Ae = 0.224 (top) and g/Ae = 
0.44 (bottom). In each case, curves corresponding to even and 
odd particle number are shown, the latter being displayed 
with additional filled circles. In the insets, relative error (in 
percent) on the total energy with respect to the exact solution 
made in the BCS and PBCS functionals is shown. 



A increases (ii) The additional approximations made to 
obtain the functional that lead to an increase of the devi- 
ation compared to PBCS as A — > 0. Nevertheless, we see 
from this comparison that the PBCS based functional is 
much more competitive than the BCS theory and is ex- 
pected to be much easier to implement than PBCS itself. 



B. Application to randomly spaced levels 

As a final illustration of the functional theory appli- 
cation, we consider here a set of randomly spaced levels. 
Following ref. (36l . l37j , an ensemble of random spectrum 
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is generated by the central eigenvalues of a 2A x 2A ran- 
dom matrix. Thus, the set of energy levels belongs to 
the Gaussian Orthogonal Ensemble, see ref. [38|. The 
renormalization proposed in ref. [37j | is performed where 

e -> 1/2tt 4Asin. _1 (e/ViAj - e\J \A - e 2 (41) 

so that the average level energy spacing is of the order of 
unity. As an illustration, evolution of average condensa- 
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FIG. 9. (Color online) Evolution of the average condensation 
energy and its statistical fluctuation (displayed by errorbars) 
as a function of g/Ae for the the PBCS based functional (blue 
filled circles), the BCS functional (green dot line) for A — 41 
(top) A = 16 (bottom) . The exact solution (red solid line) for 
an equidistant level spacing Aei of unit is shown as reference. 

tion energy and its statistical fluctuation as a function of 
g/Ae are shown in figure[§]obtained with the PBCS func- 
tional (filled circle) and the BCS functional (dot line). 
Again and as expected, the new functional matches the 
reference result of the exact solution for an equidistant 
level spacing Ae = 1. This last application illustrates 
that the method can be applied to systems with various 
level-densities. 



IV. CONCLUSION 

In this work, quasi-particle states projected onto good 
particle number are used as a starting point to propose 
new functionals dedicated to pairing correlations. The 
properties of projected states are first reviewed. These 



properties are then used to get a functional of occu- 
pation numbers and natural orbitals of the trial wave- 
function. The new functional is benchmarked with the 
pairing Hamiltonian either with equidistant or with ran- 
domly distributed single-particle energies for even and 
odd systems. In all cases, a very good agreement with 
the exact result is obtained showing great improvement 
compared to the BCS theory. Origins of the remaining 
deviations are discussed. 

The possibility to use a new functional accounting for 
particle number conservation opens new perspectives for 
the study of mesoscopic systems where pairing plays an 
important role. One may for instance anticipate new ap- 
plication for thermodynamics or dynamics where direct 
projection are too complex to provide a practical tool. In 
addition, this might also be a tool of choice avoiding re- 
cent difficulties encountered in nuclear structure studies 
(see for instance [HI). 



Appendix A: Further properties of Ik, Ik(i), ■•■ 

In this appendix, properties of the overlaps defined in 
Eqs. (TTTj) are further developed. The discussion below 
is especially useful to make connection with recent works 
and between PBCS and BCS states. Using expression 
(jTTJ) for Ik and taking advantage of the recurrence rela- 
tion (fl~3"|) gives 



Ik=J2 \*i?lK-i - (N - 1) ^ N 4 /k- 2 (*). 

i i 

In a similar way, lK-i{i) can be expressed in terms of 
Ik -2 through 

I K -2(i) = Ik-2 -{K- 2)\xi\ 2 I K -z{i), 
leading to 

Ik = \^\ 2i k-i ~(K-1)J2 \xi\ 4 Ik-2 + ■■■ 



Iterating this procedure K times leads to 
K 



(Al) 



where only overlaps II (with L < K ) appear in the right 
hand side and where the coefficients 

3 

are introduced. According to the above expression (|A1I) . 
any Ik can be written in a determinant form as 



Ik = 



Ai 


1 











A 2 


x x 


2 








A 3 


x 2 


Xi 


3 






Xk-1 Xk-2 Xk-3 Xk-4 ■ ■ ■ (K — 1) 
Xk Xk-i Xk-2 Xks • • ■ X\ 



(A2) 
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The same expression has been obtained by Rowe [25|, l26| 
using a completely different starting point making con- 
nection with elementary symmetric Schur polynomial. 
Besides this expression, similar to transformation be- 
tween elementary symmetric polynomials and power 
sums Xk, it is worth to mention that other relations 
linking other bases of the symmetric polynomial algebra 
exist [39( . 

The same procedure can also be followed for the differ- 
ent quantities //<:(&), leading to a form similar 
to (|A2|) where the X n have been respectively replaced by 
X n (i), X n (i,j), ... with: 



Xn(i,j)= J2 



|2n 



Using similar analogies between relations that hold in 
both PBCS and BCS case, we can also deduce: 



\N:i,j) = \N-l:i,j) = >>- 

= n 

(B3) 

Noting that the coefficient olk introduced in the text also 
verify: 



(K : i\K : i) 



(K- 1 : i\K - 1 : i) 
K : i.j l\ : i.j 
(K-l:i,j\K-l:i,j) 



(B4) 



Appendix B: Guidance from the BCS theory 

The BCS or HFB framework has played an important 
role in developing the new functional proposed in this 
work. We give here, highlights of some aspects discussed 
in the text. Let us start with a state given by Eq. ([3]). 
To connect with the PBCS notation, we write 

\N) = H(l + x k bi)\-), (Bl) 

k 

(B2) 

keeping in mind that, in the quasi- particle many-body 
case, the particle number N is only conserved in average 
and has only a meaning in the thermodynamics limit. 
In analogy with the PBCS case, we introduce the set of 
states \N — 1 : i) such that 

(N\N) = (N : i\N : i) + \xi\ 2 (N — 1 : i\N — 1 : i) 
{N\a\ ai \N) = \xi\ 2 {N-l:i\N-l:i) 
(N\b\b j \N)= X *x j (N-l:i\N-l:j). 

Starting from (|Blj) . we directly see that states verifying 
above relations also verify: 

\N:i) = \N-l:i) = --- 

=n( i +^i)i-)- 



We directly see that any of these coefficients identifies to 
1 in the BCS case. With this in mind, let us now give 
some intuition on how the BCS relation Q can eventu- 
ally be seen as a special limit of the PBCS case. Using 
different recurrence relations, it can be shown that 

n i = N\x i \*^-N(N-i)\x i \*^ 
In *n 

+ ... + (-i) Ar - 1 An N ^A 

In 

Assuming that all olk are equal to 1, gives 

n i = \x i \ 2 {l-\x i \ 2 + --- + \x i \ 2 ^}, 

which identifies with the BCS case, i.e. rij = |a;i| 2 /(l + 
|:Ej| 2 ) asf^co. 
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